// Unit Test for image class

#include "image/image.h"
#include "image/convolve.h"
#include "image/imageTransform.h"
#include "image/morphology.h"

#include "unitTest.h"

#include "fastTiming.h"

using namespace simple;

namespace simple_test{

    //
    // Test Functions
    //

    template<class T1>
    bool testImage(const Image<T1>& image, T1 value)
    {
        typename Image<T1>::row_iterator i;
        long w = image.width ();
        const T1* p;
        for (i=image.row_begin(); i != image.row_end(); i++) {
            p = i->p;
            for (int x=0; x<w; x++)
                if (*p++ != value)
                    return false;
        }

        return true;
    }


    template<class T>
    void dumpImage(const Image<T>& image)
    {
        typename Image<T>::row_iterator i;
        long w = image.width ();
        const T* p;
        for (i=image.row_begin(); i != image.row_end(); i++) {
            p = i->p;
            for (int x=0; x<w; x++)
                std::cout << (int) *p++ << " ";
            std::cout << std::endl;
        }
    }

    UTFUNC (templatePel8)
    {
        setDescription ("Template Pel8 instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<Pel8> image1 (boundary);
        Image<Pel8> image2 (boundary);
        image1.set (1);
        image2.set (2);
        VERIFY (testImage(image1, (Pel8)1));
        VERIFY (testImage(image2, (Pel8)2));

        Image<Pel8> image3;

        image1.mul (3);
        VERIFY (testImage(image1, (Pel8)3));
        image1.div (3);
        VERIFY (testImage(image1, (Pel8)1));
        image1.add (2);
        VERIFY (testImage(image1, (Pel8)3));
        image1.sub (1);
        VERIFY (testImage(image1, (Pel8)2));

        image3 = image1 - image2;
        VERIFY (testImage(image3, (Pel8)0));

        image3 = image1 + image2;
        VERIFY (testImage(image3, (Pel8)4));

        image2 = unsharpMask<Pel32s> (image1, 2.0);
        laplacian3x3<Pel32s> (image1, image2);
        lopass3x3<Pel32s> (image1, image2);
        hipass3x3<Pel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<Pel32s> (image1, kernel, 3, 9, image2);
    }

    UTFUNC (templatePel82)
    {
        setDescription ("Template Pel82 instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<Pel8> image1 (boundary);
        Image<Pel8> image2 (boundary);
        Image<Pel8> image3;

        image1.set (1);
        image1.mul (1);
        image1.div (1);
        image1.add (1);
        image1.sub (1);

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<Pel8> (image1, 2.0);
        laplacian3x3<Pel8> (image1, image2);
        lopass3x3<Pel8> (image1, image2);
        hipass3x3<Pel8> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<Pel8> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templatePel16)
    {
        setDescription ("Template Pel16 instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<Pel16> image1 (boundary);
        Image<Pel16> image2 (boundary);
        Image<Pel16> image3;

        image1.set (1);
        image1.mul (1);
        image1.div (1);
        image1.add (1);
        image1.sub (1);

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<Pel32s> (image1, 20);
        laplacian3x3<Pel32s> (image1, image2);
        lopass3x3<Pel32s> (image1, image2);
        hipass3x3<Pel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<Pel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templatePel32s)
    {
        setDescription ("Template Pel32s instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<Pel32s> image1 (boundary);
        Image<Pel32s> image2 (boundary);
        Image<Pel32s> image3;

        image1.set (1);
        image1.mul (1);
        image1.div (1);
        image1.add (1);
        image1.sub (1);

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<Pel32s> (image1, 2.0);
        laplacian3x3<Pel32s> (image1, image2);
        lopass3x3<Pel32s> (image1, image2);
        hipass3x3<Pel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<Pel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templateRGB)
    {
        setDescription ("Template RGB instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<RGB> image1 (boundary);
        Image<RGB> image2 (boundary);
        Image<RGB> image3;

        image1.set (RGB(1));
        image1.mul (RGB(1));
        image1.div (RGB(1));
        image1.add (RGB(1));
        image1.sub (RGB(1));

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<RGBPel32s> (image1, 2.0);
        laplacian3x3<RGBPel32s> (image1, image2);
        lopass3x3<RGBPel32s> (image1, image2);
        hipass3x3<RGBPel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<RGBPel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templateRGB2)
    {
        setDescription ("Template RGB2 instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<RGB> image1 (boundary);
        Image<RGB> image2 (boundary);
        Image<RGB> image3;

        image1.set (RGB(1));
        image1.mul (RGB(1));
        image1.div (RGB(1));
        image1.add (RGB(1));
        image1.sub (RGB(1));

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<RGB> (image1, 2.0);
        laplacian3x3<RGB> (image1, image2);
        lopass3x3<RGB> (image1, image2);
        hipass3x3<RGB> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<RGB> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templateRGBPel32s)
    {
        setDescription ("Template RGBPel32s instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<RGBPel32s> image1 (boundary);
        Image<RGBPel32s> image2 (boundary);
        Image<RGBPel32s> image3;

        image1.set (RGBPel32s(1));
        image1.mul (RGBPel32s(1));
        image1.div (RGBPel32s(1));
        image1.add (RGBPel32s(1));
        image1.sub (RGBPel32s(1));

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<RGBPel32s> (image1, 2.0);
        laplacian3x3<RGBPel32s> (image1, image2);
        lopass3x3<RGBPel32s> (image1, image2);
        hipass3x3<RGBPel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<RGBPel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templateClampedPel8)
    {
        setDescription ("Template ClampedPel8 instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<ClampedPel8> image1 (boundary);
        Image<ClampedPel8> image2 (boundary);
        Image<ClampedPel8> image3;

        image1.set (1);
        image1.mul (1);
        image1.div (1);
        image1.add (1);
        image1.sub (1);

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<ClampedPel32s> (image1, 2.0);
        laplacian3x3<ClampedPel32s> (image1, image2);
        lopass3x3<ClampedPel32s> (image1, image2);
        hipass3x3<ClampedPel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<ClampedPel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC (templateClampedPel32s)
    {
        setDescription ("Template ClampedPel32s instantiation");

        Rect boundary (0, 0, 16, 16);

        Image<ClampedPel32s> image1 (boundary);
        Image<ClampedPel32s> image2 (boundary);
        Image<ClampedPel32s> image3;

        image1.set (1);
        image1.mul (1);
        image1.div (1);
        image1.add (1);
        image1.sub (1);

        image3 = image1 - image2;
        image3 = image1 + image2;

        image2 = unsharpMask<ClampedPel32s> (image1, 2.0);
        laplacian3x3<ClampedPel32s> (image1, image2);
        lopass3x3<ClampedPel32s> (image1, image2);
        hipass3x3<ClampedPel32s> (image1, image2);

        char kernel[] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
        convolve<ClampedPel32s> (image1, kernel, 3, 9, image2);
    }


    UTFUNC(RGB)
    {
        setDescription ("RGB tests");

        Pel8 b = 250;
        ClampedPel8 bc = 250;

        b += 10;
        bc += 10;

        std::cout << "b = " << (long)b << std::endl;
        std::cout << "bc = " << bc << std::endl;

        RGB rgb (250);
        rgb += 10;
        std::cout << "rgb = " << rgb << std::endl;

        RGBTmpl<Pel8> rgb1 (1,2,3);
        std::cout << rgb1 << std::endl;
        rgb1 += rgb1;
        std::cout << rgb1 << std::endl;

        RGBTmpl<Pel32s> rgb2 (4,5,6);
        std::cout << rgb2 << std::endl;

        rgb2 += rgb1;
        std::cout << rgb2 << std::endl;
    }

    /*
        {
        Rect b1 (0, 0, 16, 16);
        Image<RGB> image (b1);
        image.set (250);
        image.add (10);
        Image<RGB>::row_iterator i;
        for (i=image.row_begin(); i != image.row_end(); i++) {
            const RGB* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (int x=0; x<image.width(); x++)
                std::cout << *p++ << " ";
            std::cout << std::endl;
        }
        }
    */


    UTFUNC(ctor)
    {
        setDescription ("Constructor and simple accessor tests");

        Rect boundary (0, 0, 640, 480);

        Image<Pel8> Pel8Image1;
        Image<Pel8> Pel8Image2 (boundary);
        Image<Pel8> Pel8Image3 (boundary, 
                                  RectImageStorage::eQuadWordAlign);

        VERIFY (Pel8Image1.isNull());
        VERIFY (!Pel8Image2.isNull());
        VERIFY (!Pel8Image3.isNull());
        VERIFY (Pel8Image2.boundary() == boundary);
        VERIFY (Pel8Image3.boundary() == boundary);
        VERIFY (Pel8Image2.bytesPerPixel() == sizeof (Pel8));
        VERIFY (Pel8Image3.bytesPerPixel() == sizeof (Pel8));
        VERIFY (Pel8Image2.xoffset() == 0);
        VERIFY (Pel8Image3.xoffset() == 0);
        VERIFY (Pel8Image2.yoffset() == 0);
        VERIFY (Pel8Image3.yoffset() == 0);



        Image<Pel32s> longImage1;
        Image<Pel32s> longImage2 (boundary);
        Image<Pel32s> longImage3 (boundary, 
                                    RectImageStorage::eQuadWordAlign);

        VERIFY (longImage1.isNull());
        VERIFY (!longImage2.isNull());
        VERIFY (!longImage3.isNull());
        VERIFY (longImage2.boundary() == boundary);
        VERIFY (longImage3.boundary() == boundary);
        VERIFY (longImage2.bytesPerPixel() == sizeof (Pel32s));
        VERIFY (longImage3.bytesPerPixel() == sizeof (Pel32s));
        VERIFY (longImage2.xoffset() == 0);
        VERIFY (longImage3.xoffset() == 0);
        VERIFY (longImage2.yoffset() == 0);
        VERIFY (longImage3.yoffset() == 0);



        Image<RGB> RGBImage1;
        Image<RGB> RGBImage2 (boundary);
        Image<RGB> RGBImage3 (boundary, 
                                    RectImageStorage::eQuadWordAlign);

        VERIFY (RGBImage1.isNull());
        VERIFY (!RGBImage2.isNull());
        VERIFY (!RGBImage3.isNull());
        VERIFY (RGBImage2.boundary() == boundary);
        VERIFY (RGBImage3.boundary() == boundary);
        VERIFY (RGBImage2.bytesPerPixel() == sizeof (RGB));
        VERIFY (RGBImage3.bytesPerPixel() == sizeof (RGB));
        VERIFY (RGBImage2.xoffset() == 0);
        VERIFY (RGBImage3.xoffset() == 0);
        VERIFY (RGBImage2.yoffset() == 0);
        VERIFY (RGBImage3.yoffset() == 0);
    }


    UTFUNC(setPel8)
    {
        setDescription ("set Pel8");

        int i;
        Rect boundary (0, 0, 1024, 1024);

        Image<Pel8> Pel8Image1 (boundary);
        for (i=0; i<100; i++)
            Pel8Image1.set (0);

        Image<Pel8> Pel8Image2 = Pel8Image1;
        VERIFY (Pel8Image1 == Pel8Image2);
    }


    UTFUNC(setlong)
    {
        setDescription ("set Pel32s");

        int i;
        Rect boundary (0, 0, 1024, 1024);

        Image<Pel32s> longImage1 (boundary);
        for (i=0; i<100; i++)
            longImage1.set (0);

        Image<Pel32s> longImage2 = longImage1;
        VERIFY (longImage1 == longImage2);
    }

    UTFUNC(setRGB)
    {
        setDescription ("set RGB");

        int i;
        Rect boundary (0, 0, 1024, 1024);

        Image<RGB> RGBImage1 (boundary);
        for (i=0; i<100; i++)
            RGBImage1.set (RGB(0));

        Image<RGB> RGBImage2 = RGBImage1;
        VERIFY (RGBImage1 == RGBImage2);
    }


    UTFUNC(mixed)
    {
        setDescription ("Mixed");

        Rect boundary (0, 0, 1024, 1024);

        Image<Pel8> Pel8Image1 (boundary);
        Pel8Image1.set (1);

        Image<Pel32s> longImage1 (boundary);
        longImage1.set (0);

        add (Pel8Image1, longImage1);
        VERIFY (testImage(longImage1, (Pel32s)1));

        longImage1 += Pel8Image1;
        VERIFY (testImage(longImage1, (Pel32s)2));

        sub (Pel8Image1, longImage1);
        VERIFY (testImage(longImage1, (Pel32s)1));

        longImage1 -= Pel8Image1;
        VERIFY (testImage(longImage1, (Pel32s)0));

        Image<Pel32s> longImage2 (boundary);
        longImage1.set (1);
        longImage2.set (1);
        VERIFY (testImage(longImage2, (Pel32s)1));

        // Copy unlike images
        Image<Pel8> Pel8Image2;
        copy (longImage2, Pel8Image2);

        VERIFY (longImage1 == longImage2);
        VERIFY (Pel8Image1 == Pel8Image2);
        VERIFY (testImage(Pel8Image2, (Pel8)1));

        // Copy like images
        Image<Pel8> Pel8Image3;
        copy (Pel8Image2, Pel8Image3);
        VERIFY (Pel8Image2 == Pel8Image3);
        VERIFY (Pel8Image2.base() != Pel8Image3.base());
        VERIFY (testImage(Pel8Image3, (Pel8)1));

        // Simple copy
        Image<Pel8> Pel8Image4 = copy (Pel8Image3);
        VERIFY (Pel8Image3 == Pel8Image4);
        VERIFY (Pel8Image3.base() != Pel8Image4.base());

        // Boundary behavior (wrapping)
        Pel8Image1.set (255);
        Pel8Image2.set (1);
        VERIFY (testImage(Pel8Image1, (Pel8)255));
        Pel8Image1 += Pel8Image2;
        VERIFY (testImage(Pel8Image1, (Pel8)0));

        Pel8Image1.set (0);
        Pel8Image1 -= Pel8Image2;
        VERIFY (testImage(Pel8Image1, (Pel8)255));

    }


    UTFUNC(clampedMath)
    {
        setDescription ("Mixed");

        Rect boundary (0, 0, 1024, 1024);

        Image<ClampedPel8> Pel8Image1 (boundary);
        Image<ClampedPel8> Pel8Image2 (boundary);
        Pel8Image1.set (254);
        Pel8Image2.set (1);
        
        VERIFY (testImage(Pel8Image1, (ClampedPel8)254));
        VERIFY (testImage(Pel8Image2, (ClampedPel8)1));

        Pel8Image1 += Pel8Image2;
        VERIFY (testImage(Pel8Image1, (ClampedPel8)255));

        // Should clamp at 255
        Pel8Image1 += Pel8Image2;
        VERIFY (testImage(Pel8Image1, (ClampedPel8)255));

        Pel8Image1 -= Pel8Image2;
        VERIFY (testImage(Pel8Image1, (ClampedPel8)254));
    }


    UTFUNC(mixedrgb)
    {
        setDescription ("Mixed rgb");

        Rect boundary (0, 0, 1024, 1024);

        Image<RGB> image1 (boundary);
        RGB rgb (10, 20, 30);
        image1.set (rgb);

        Image<Pel8> image2 (boundary);
        copy (image1, image2);

        // Verify each pixel
        VERIFY (true);
        Image<Pel8>::row_iterator i;
        long w = image2.width ();
        const Pel8* p;
        for (i=image2.row_begin(); i != image2.row_end(); i++) {
            p = i->p;
            for (int x=0; x<w; x++)
                if (*p++ != 20)
                    VERIFY (!"Pixel verification error 1");
        }

        image1.set (rgb);
        image2.set (5);
        image1 += image2;
        RGB rgb2 (15, 25, 35);
        VERIFY (testImage (image1, rgb2));

    }


    UTFUNC(align)
    {
        setDescription ("align");

        // This image is NOT aligned
        Rect boundary (0, 0, 30, 30);
        Image<Pel8> image1 (boundary);

        Image<Pel8> align0 = image1.align (RectImageStorage::eNoAlign);
        Image<Pel8> align2 = image1.align (RectImageStorage::eWordAlign);
        Image<Pel8> align4 = image1.align (RectImageStorage::eDoubleWordAlign);
        Image<Pel8> align8 = image1.align (RectImageStorage::eQuadWordAlign);

        VERIFY (image1.rowAddress(0) == align0.rowAddress(0));
        VERIFY (align0.alignment() == RectImageStorage::eNoAlign);
        VERIFY (align2.bestAlignment() == RectImageStorage::eDoubleWordAlign ||
                        align2.bestAlignment() == RectImageStorage::eQuadWordAlign ||
                        align2.bestAlignment() == RectImageStorage::eWordAlign);
        VERIFY (align4.bestAlignment() == RectImageStorage::eDoubleWordAlign ||
                        align4.bestAlignment() == RectImageStorage::eQuadWordAlign);
        VERIFY (align8.bestAlignment() == RectImageStorage::eQuadWordAlign);

        // Do the same thing using our duplicate function
        align0 = duplicate (image1, RectImageStorage::eNoAlign);
        align2 = duplicate (image1, RectImageStorage::eWordAlign);
        align4 = duplicate (image1, RectImageStorage::eDoubleWordAlign);
        align8 = duplicate (image1, RectImageStorage::eQuadWordAlign);

        VERIFY (image1.rowAddress(0)    != align0.rowAddress(0));
        VERIFY (align0.alignment()      == RectImageStorage::eNoAlign);
        VERIFY (align2.bestAlignment()  == RectImageStorage::eDoubleWordAlign ||
                align2.bestAlignment()  == RectImageStorage::eQuadWordAlign ||
                align2.bestAlignment()  == RectImageStorage::eWordAlign);
        VERIFY (align4.bestAlignment()  == RectImageStorage::eDoubleWordAlign ||
                align4.bestAlignment()  == RectImageStorage::eQuadWordAlign);
        VERIFY (align8.bestAlignment()  == RectImageStorage::eQuadWordAlign);

    }

    UTFUNC(trim)
    {
        setDescription ("trim");

        Rect boundary (0, 0, 100, 100);
        Image<Pel8> image1 (boundary);

        Rect window (10, 20, 30, 40);
        Image<Pel8> image2 = image1;
        image2.window (window);

        VERIFY (image2.x0()     == 10);
        VERIFY (image2.y0()     == 20);
        VERIFY (image2.width()  == 30);
        VERIFY (image2.height() == 40);
        VERIFY (image2.ref()    == 2);

        VERIFY (image2.storage().xoffset() == 10);
        VERIFY (image2.storage().yoffset() == 20);

        image2.trim ();

        VERIFY (image2.x0()         == 10);
        VERIFY (image2.y0()         == 20);
        VERIFY (image2.width()    == 30);
        VERIFY (image2.height() == 40);
        VERIFY (image2.ref()        == 1);

        VERIFY (image2.storage().xoffset() == 0);
        VERIFY (image2.storage().yoffset() == 0);

    }


    UTFUNC(rebase)
    {
        setDescription ("rebase");

        Rect boundary (0, 0, 100, 100);
        Image<Pel8> image1 (boundary);

        Rect window (10, 20, 30, 40);
        Image<Pel8> image2 = image1;
        image2.window (window);

        VERIFY (image2.x0()         == 10);
        VERIFY (image2.y0()         == 20);
        VERIFY (image2.width()    == 30);
        VERIFY (image2.height() == 40);
        VERIFY (image2.ref()        == 2);

        VERIFY (image2.storage().xoffset() == 10);
        VERIFY (image2.storage().yoffset() == 20);

        image2.rebase ();

        VERIFY (image2.x0()         == 0);
        VERIFY (image2.y0()         == 0);
        VERIFY (image2.width()    == 30);
        VERIFY (image2.height() == 40);
        VERIFY (image2.ref()        == 1);

        VERIFY (image2.storage().xoffset() == 0);
        VERIFY (image2.storage().yoffset() == 0);

    }


    UTFUNC(add)
    {
        setDescription ("add");

        Rect boundary1 (0, 0, 640, 480);
        Rect boundary2 (1, 2, 638, 476);
        Rect boundary3 (2, 4, 636, 472);

        Image<Pel8> Pel8Image1 (boundary1);
        Pel8Image1.set (1);

        Image<Pel8> Pel8Image2 (boundary2);
        Pel8Image2.set (2);

        Image<Pel8> Pel8Image3;

        // No destination specified
        add (Pel8Image1, Pel8Image2, Pel8Image3);

        Rect overlap = boundary1.intersect (boundary2);
        VERIFY (Pel8Image3.boundary() == overlap);
        VERIFY (testImage (Pel8Image3, (Pel8)3));


        // Destination specified
        Pel8Image3 = Image<Pel8> (boundary3);
        Pel8Image3.set (100);

        add (Pel8Image1, Pel8Image2, Pel8Image3);
        VERIFY (Pel8Image3.boundary() == boundary3);
        VERIFY (testImage (Pel8Image3, (Pel8)3));


        // operator+ syntax. Additional causes wrapping
        Pel8Image1.set (255);
        Pel8Image2.set (2);
        Pel8Image3 = Pel8Image1 + Pel8Image2;

        VERIFY (Pel8Image3.boundary() == overlap);
        VERIFY (testImage (Pel8Image3, (Pel8)1));

        // operator+ syntax. Additional causes clamping
        Image<ClampedPel8> Pel8Image1c (boundary1);
        Image<ClampedPel8> Pel8Image2c (boundary2);
        Image<ClampedPel8> Pel8Image3c;

        Pel8Image1c.set (254);
        Pel8Image2c.set (2);
        Pel8Image3c = Pel8Image1c + Pel8Image2c;

        VERIFY (Pel8Image3c.boundary() == overlap);
        VERIFY (testImage (Pel8Image3c, (ClampedPel8)255));

    }


    UTFUNC (convolve)
    {
        setDescription ("convolve");

        Rect boundary1 (0, 0, 640, 480);

        Image<Pel8> Pel8Image1 (boundary1);
        Image<Pel8> Pel8Image2 (boundary1);
        Pel8Image1.set (1);
        Pel8Image2.set (0);

        char kernel[] = { (char)-1, (char)-1, (char)-1, (char)-1, 8, (char)-1, (char)-1, (char)-1, (char)-1 };

        convolve<Pel32s> (Pel8Image1, kernel, 3, 1, Pel8Image2);

        Image<Pel8> Pel8Image3 (boundary1);
        Pel8Image3.set (0);
        laplacian3x3<Pel32s> (Pel8Image1, Pel8Image3);
        std::cout << "Pel8Image1 = " << (long) Pel8Image1.getPixel(0,0) << std::endl;
        std::cout << "Pel8Image2 = " << (long) Pel8Image2.getPixel(0,0) << std::endl;
        std::cout << "Pel8Image3 = " << (long) Pel8Image3.getPixel(0,0) << std::endl;

        VERIFY (Pel8Image2 == Pel8Image3);

        Image<Pel8> Pel8Image4, Pel8Image5;
        laplacian3x3<Pel32s> (Pel8Image4, Pel8Image5);
        VERIFY (Pel8Image4.isNull());
        VERIFY (Pel8Image5.isNull());

    /*
        int o;
        for (o=0; o<9; o++)
            std::cout << (int) kernel[o] << " ";
        std::cout << std::endl;

        Image<Pel8>::row_iterator i;
        for (i=Pel8Image2.row_begin(); i != Pel8Image2.row_end(); i++) {
            const Pel8* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (int x=0; x<Pel8Image2.width(); x++)
                std::cout << (Pel32s) *p++ << " ";
            std::cout << std::endl;
        }
    */

    }


    UTFUNC (laplacian3x3)
    {
        setDescription ("laplacian3x3");

        Rect boundary1 (0, 0, 640, 480);

        Image<Pel8> Pel8Image1 (boundary1);
        Image<Pel8> Pel8Image2 (boundary1);
        Pel8Image1.set (1);

        laplacian3x3<Pel32s> (Pel8Image1, Pel8Image2);

        Rect boundary2 (0, 0, 638, 478);
        Image<Pel8> Pel8Image3 (boundary2);
        Pel8Image3.set (0);

        Image<Pel8> Pel8Image4 = Pel8Image2;
        Pel8Image4.window (boundary2);
        VERIFY (Pel8Image3 == Pel8Image4);
        
    /*
        int o;
        for (o=0; o<9; o++)
            std::cout << (int) kernel[o] << " ";
        std::cout << std::endl;

        Image<Pel8>::row_iterator i;
        for (i=Pel8Image2.row_begin(); i != Pel8Image2.row_end(); i++) {
            const Pel8* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (int x=0; x<Pel8Image2.width(); x++)
                std::cout << (Pel32s) *p++ << " ";
            std::cout << std::endl;
        }
    */

    }


    UTFUNC (laplacian3x3dump)
    {
        setDescription ("laplacian3x3 dump");

        Rect boundary1 (0, 0, 16, 16);

        Image<Pel8> Pel8Image1 (boundary1);
        Image<Pel8> Pel8Image2 (boundary1);
        Pel8Image1.set (200);
        Pel8Image2.set (2);
        Pel8Image1.setPixel (Point(8,8), 2);

        laplacian3x3<long> (Pel8Image1, Pel8Image2);
        VERIFY (Pel8Image2.boundary() == boundary1);

        Image<Pel8> Pel8Image3;
        laplacian3x3<long> (Pel8Image1, Pel8Image3);

        Rect boundary2 = boundary1;
        boundary2.expand (-1, -1);
        VERIFY (Pel8Image3.boundary() == boundary2);

        Image<RGB> rgbImage1 (boundary1);
        Image<RGB> rgbImage2;
        rgbImage1.set (RGB(200));
        rgbImage1.setPixel (Point(8,8), RGB(2));

        laplacian3x3<RGBPel32s> (rgbImage1, rgbImage2);
        VERIFY (rgbImage2.boundary() == boundary2);


        Image<Pel8>::row_iterator i;
        for (i=Pel8Image2.row_begin(); i != Pel8Image2.row_end(); i++) {
            const Pel8* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (unsigned int x=0; x<Pel8Image2.width(); x++)
                std::cout << (long) *p++ << " ";
            std::cout << std::endl;
        }

        for (i=Pel8Image3.row_begin(); i != Pel8Image3.row_end(); i++) {
            const Pel8* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (unsigned int x=0; x<Pel8Image3.width(); x++)
                std::cout << (long) *p++ << " ";
            std::cout << std::endl;
        }

        Image<RGB>::row_iterator i1;
        for (i1=rgbImage2.row_begin(); i1 != rgbImage2.row_end(); i1++) {
            const RGB* p = i1->p;
            std::cout << i1->x << ", " << i1->y << ": ";
            for (unsigned int x=0; x<rgbImage2.width(); x++)
                std::cout << *p++ << " ";
            std::cout << std::endl;
        }


    }

    UTFUNC (laplacian3x3RGB)
    {
        setDescription ("laplacian3x3 RGB");

        Rect boundary1 (0, 0, 640, 480);

        Image<RGB> RGBImage1 (boundary1);
        Image<RGB> RGBImage2 (boundary1);
        RGBImage1.set (RGB(1));

        laplacian3x3<RGBPel32s> (RGBImage1, RGBImage2);

        Rect boundary2 (0, 0, 638, 478);
        Image<RGB> RGBImage3 (boundary2);
        RGBImage3.set (RGB(0));

        Image<RGB> RGBImage4 = RGBImage2;
        RGBImage4.window (boundary2);
        VERIFY (RGBImage3 == RGBImage4);
    }



    UTFUNC (laplacian3x3double)
    {
        setDescription ("laplacian3x3 double");

        Rect boundary1 (0, 0, 640, 480);

        Image<double> doubleImage1 (boundary1);
        Image<double> doubleImage2 (boundary1);
        doubleImage1.set (1.);

        laplacian3x3<double> (doubleImage1, doubleImage2);

        Rect boundary2 (0, 0, 638, 478);
        Image<double> doubleImage3 (boundary2);
        doubleImage3.set (0);

        Image<double> doubleImage4 = doubleImage2;
        doubleImage4.window (boundary2);
        VERIFY (doubleImage3 == doubleImage4);
    }


    UTFUNC (unsharpmask)
    {
        setDescription ("unsharp mask");

        Rect boundary1 (0, 0, 640, 480);

        Image<RGB> image1 (boundary1);
        Image<RGB> image2 (boundary1);
        image1.set (RGB(1));

        image2 = unsharpMask<RGBPel32s> (image1, 2.0);

    }


    UTFUNC (thumbnail)
    {
        setDescription ("thumbnail");

        Rect boundary1 (0, 0, 8, 8);
        Image<RGB> image (boundary1);

        Image<RGB>::row_iterator i;
        Pel8 count = 0;
        for (i=image.row_begin(); i != image.row_end(); i++) {
            RGB* p = i->p;
            for (unsigned int x=0; x<image.width(); x++)
                *p++ = RGB(count++);
        }

        std::cout << "Original:" << std::endl;
        for (i=image.row_begin(); i != image.row_end(); i++) {
            const RGB* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (unsigned int x=0; x<image.width(); x++)
                std::cout << *p++ << " ";
            std::cout << std::endl;
        }

        Image<RGB> tnail = thumbnail<RGBPel32s> (image, 2);

        std::cout << "Thumbnail:" << std::endl;
        for (i=tnail.row_begin(); i != tnail.row_end(); i++) {
            const RGB* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (unsigned int x=0; x<tnail.width(); x++)
                std::cout << *p++ << " ";
            std::cout << std::endl;
        }

        // This will cause wrapping since intermediate quantity isn't big enough
        Image<RGB> tnail2 = thumbnail<RGB> (image, 3);

        std::cout << "Thumbnail Pel8 average:" << std::endl;
        for (i=tnail2.row_begin(); i != tnail2.row_end(); i++) {
            const RGB* p = i->p;
            std::cout << i->x << ", " << i->y << ": ";
            for (unsigned int x=0; x<tnail2.width(); x++)
                std::cout << *p++ << " ";
            std::cout << std::endl;
        }

    }


    UTFUNC (morphology)
    {
        setDescription ("morphology");

        Rect boundary1 (0, 0, 8, 8);
        Image<RGB> image (boundary1);

        RGB black (0);
        RGB white (255);

        // Create a cross image
        image.set (black);
        image.setPixel (4, 4, white);
        image.setPixel (3, 4, white);
        image.setPixel (5, 4, white);
        image.setPixel (4, 3, white);
        image.setPixel (4, 5, white);

        Image<RGB> dst;
        erodeCross (image, dst);

        // Verify it again the correct output
        Rect boundary2 (1, 1, 6, 6);
        Image<RGB> correct (boundary2);
        correct.set (black);
        correct.setPixel (4, 4, white);
        VERIFY (dst == correct);
    }


    UTFUNC (timing)
    {
        int i;

        setDescription ("timing");

        // copy performance compared to stl version
        Rect boundary (0, 0, 1024, 1024);
        Image<Pel8> Pel8Image1 (boundary);
        Pel8Image1.set (1);

        HiResElapsedTime t1;
        Image<Pel8> Pel8Image2;
        for (i=0; i<10; i++)
            Pel8Image2 = copy (Pel8Image1);
        double t1msec = t1.msec () / 10.;

        HiResElapsedTime t2;
        Image<Pel8> Pel8Image3;
        for (i=0; i<10; i++)
            Pel8Image3 = copy_stl (Pel8Image1);
        double t2msec = t2.msec () / 10.;

        VERIFY (Pel8Image1 == Pel8Image2);
        VERIFY (Pel8Image1 == Pel8Image3);

        std::cout << "1024x1024 copy in " << t1msec << " msec" << std::endl;
        std::cout << "1024x1024 copy_stl in " << t2msec << " msec" << std::endl;

        // Test generic convolution vs. faster version
        char kernel[] = {-1, -1, -1, -1, 8, -1, -1, -1, -1};
        Pel8Image1.set (1);

        HiResElapsedTime t3;
        Image<Pel8> laplace1;
        for (i=0; i<10; i++) {
            laplace1 = Image<Pel8> (boundary);
            laplace1.set (0);
            convolve<Pel32s> (Pel8Image1, kernel, 3, 1, laplace1);
        }
        double t3msec = t3.msec () / 10.;
     
        HiResElapsedTime t4;
        Image<Pel8> laplace2;
        for (i=0; i<10; i++) {
         laplace2 = Image<Pel8> (boundary);
         laplace2.set (0);
         laplacian3x3<Pel32s> (Pel8Image1, laplace2);
        }
        double t4msec = t4.msec () / 10.;

        std::cout << "1024x1024 pel8 laplace via convolve() in " << t3msec << " msec" << std::endl;
        std::cout << "1024x1024 pel8 laplace in " << t4msec << " msec" << std::endl;

        Image<RGB>rgbImage1 (boundary);
        rgbImage1.set (RGB(1,1,1));

        HiResElapsedTime t5;
        Image<Pel8> laplace3;
        for (i=0; i<10; i++) {
            laplace3 = Image<RGB> (boundary);
            rgbImage1.set (RGB(0,0,0));
            convolve<RGBPel32s> (rgbImage1, kernel, 3, 1, laplace3);
        }
        double t5msec = t5.msec () / 10.;
     
        HiResElapsedTime t6;
        Image<Pel8> laplace4;
        for (i=0; i<10; i++) {
            laplace4 = Image<RGB> (boundary);
            rgbImage1.set (RGB(0,0,0));
            laplacian3x3<RGBPel32s> (rgbImage1, laplace4);
        }
        double t6msec = t6.msec () / 10.;

        std::cout << "1024x1024 rgb laplace via convolve() in " << t5msec << " msec" << std::endl;
        std::cout << "1024x1024 rgb laplace in " << t6msec << " msec" << std::endl;

        VERIFY (laplace1 == laplace2);
        VERIFY (laplace3 == laplace4);
    }

}

